#name : updateRsSub.pm
#author : Alexandre DUVAL
#mail : alexandre.duval@nantes.fr
#create date : 07.21.2010
#modified : 11.04.2010	

#This library contains the Principal subroutines of updateRs.pl
  
use strict;
use DBI;

package updateRsSub;

use updteRsUtils qw(:DEFAULT);
use Exporter;
use vars qw(@ISA @EXPORT);
@ISA=qw(Exporter);
@EXPORT=qw(LOAD_DB LOAD_RsMERGE COMPARE1000G STRAND UPDATERSNAME CONTROL_BUILD CONTROL_DUPLICATE_1000G);

###################################################################################################################################################################
sub LOAD_DB
{
	my ($HashRef,$dir,$file1,$file2)=@_;
	
	my %Hashtmp=();
	my $fileDB1=$file1;
	open(DB1, $fileDB1) or die "Cannot open database file $fileDB1 : $!\n";
	while (my $lineDB1=<DB1>)
	{
		chomp ($lineDB1);
		my @line=split(/ /,$lineDB1); #chrom pos combinated_rate GeneticMap Bloc
		$Hashtmp{$line[0],$line[1]}[0]=$line[3];
		$Hashtmp{$line[0],$line[1]}[1]=$line[4];
		
	}


	my $fileDB2=$dir."/".$file2;
	open(DB2, $fileDB2) or die "Cannot open database file $fileDB2 : $!\n";
	
	while (my $lineDB2= <DB2>) 
	{
		chomp ($lineDB2);
		my @row=split(/ /,$lineDB2); #$chr rsId posgen a1 a2 freq_a1
		$$HashRef{$row[1]}[0]=1; #Exist ?
		$$HashRef{$row[1]}[1]=$row[0]; #chr
		$$HashRef{$row[1]}[2]=$row[1]; #rsid
		$$HashRef{$row[1]}[3]=$row[2]; #pos
		$$HashRef{$row[1]}[4]=$row[3]; #a1
		$$HashRef{$row[1]}[5]=$row[4]; #a2
		
		$$HashRef{$row[1]}[6]{$row[3]}=$row[5]; #freq a1
		my $freqa2=1-$row[5];
		$$HashRef{$row[1]}[6]{$row[4]}=sprintf("%1.3f",$freqa2); #freq a2

		if ( exists($Hashtmp{$row[0],$row[2]}[0]) )	
		{ $$HashRef{$row[1]}[7]=$Hashtmp{$row[0],$row[2]}[0] } #genetic_position
		else { $$HashRef{$row[1]}[7]="NA" }
		if ( exists($Hashtmp{$row[0],$row[2]}[1]) )	
		{ $$HashRef{$row[1]}[8]=$Hashtmp{$row[0],$row[2]}[1] } #block name
		else { $$HashRef{$row[1]}[8]="NA" }
		
	
	}
	close(DB); 

}
########################################################################################
sub FREQ_DATA
{
	my($file,$hashRef)=@_;
	open(FREQ,"$file") or die "cannot open FREQDATA file : $!\n";
	

	while (my $line=<FREQ>)
	{	
		$line=~s/ +/ /g;
		my @value=split(/ /,$line); #	CHR         SNP   A1   A2          MAF  NCHROBS (nota : the file start with space so the first index of the table value  is 1) 

		$$hashRef{$value[2],$value[3]}=$value[5]; #freq A1
		my $freqA2=1-$value[5];
		$$hashRef{$value[2],$value[4]}=sprintf("%1.3f",$freqA2); #freq A2
	}
	
}
#########################################################################################
sub UPDATE_RM
{
	my ($path,$pre,$post,$threeshold,$outfile)=@_;
	open(OUT,">",$outfile) or die "Cannot open output file:$!\n";

	for (my $i=1;$i<23;$i++)
	{
		my $file=$path."/".$pre."chr".$i.$post;
		open(FILE,$file) or die "Cannot open input file $file :$!\n";
		my $maxRate=-1;
		my $max=0;
		my $bloc=1;
		while (my $ligne=<FILE>)
		{
			chomp($ligne);
			my($position,$rate,$genMap)=split(/ /,$ligne);
			if ($position=="position")
			{	
				print OUT "chromosome".$ligne." block_partition\n";
				next;
			}
	
			if (($rate >= $threeshold))
			{
				if($rate > $maxRate) { $maxRate=$rate;}	
				print OUT $i." ".$position." ".$rate." ".$genMap." ".$i."-HS\n";
			}
			elsif (($maxRate > $rate))
			{
				$maxRate=-1;
				$bloc++;
	
				if($rate < $threeshold)
				{
					print OUT $i." ".$position." ".$rate." ".$genMap." ".$i."-".$bloc."\n";
				}
				else{print OUT $i." ".$position." ".$rate." ".$genMap." ".$i."-HS\n";}
			}
			else {print OUT $i." ".$position." ".$rate." ".$genMap." ".$i."-".$bloc."\n";}
	
		 }

	close(FILE);
	}
	close(OUT);

}

###################################################################################################################################################################
sub LOAD_RsMERGE
{
        my ($RsMergeFile, $newNameRef, $existMergeRef)=@_;

	OPEN ($RsMergeFile, \*FHI);
	#open (FHI, $RsMergeFile) or die "Cannot open $RsMergeFile : $!\n";  
	while (my $lineMerge = <FHI>) 
	{
	    	$lineMerge =~ s/\012//g;
	    	$lineMerge =~ s/\015//g;
	    	my ($rsHigh,$rsLow,$build_id,$orien,$create_time,$last_updated_time,$rsCurrent,$orien2Current) = split(/\t/,$lineMerge);   
	    	my $oldName="rs".$rsHigh;
	    	$$newNameRef{$oldName}= "rs".$rsLow; 
	    	$$existMergeRef{$oldName}=1;	     	         
	}
 
    	close(FHI);
}
###################################################################################################################################################################

###################################################################################################################################################################
sub MONOMORPH
{
	my($amono,$a1,$a2,$adb,$db,$rsname,$chr)=@_;

	if(($amono==1) && ($a1 eq $$db{$rsname}[4])){$a2=$$db{$rsname}[5]; $adb=$a1;}
	elsif(($amono==1) && ($a1 eq $$db{$rsname}[5])){$a2=$$db{$rsname}[4]; $adb=$a1;}
	elsif(($amono==1) && (STRAND($a1) eq $$db{$rsname}[4])){$a1=STRAND($a1); $a2=$$db{$rsname}[5]; $adb=$a1;}
	elsif(($amono==1) && (STRAND($a1) eq $$db{$rsname}[5])){$a1=STRAND($a1); $a2=$$db{$rsname}[4]; $adb=$a1;}
	elsif(($amono==2) && ($a2 eq $$db{$rsname}[4])){$a1=$$db{$rsname}[5]; $adb=$a2;}
	elsif(($amono==2) && ($a2 eq $$db{$rsname}[5])){$a1=$$db{$rsname}[4]; $adb=$a2;}
	elsif(($amono==2) && (STRAND($a2) eq $$db{$rsname}[4])){$a2=STRAND($a2); $a1=$$db{$rsname}[5]; $adb=$a2;}
	elsif(($amono==2) && (STRAND($a2) eq $$db{$rsname}[5])){$a2=STRAND($a2); $a1=$$db{$rsname}[4]; $adb=$a2;}
	else{ print EXCLUDE "$chr\t$rsname\t$a1\t$a2\t$$db{$rsname}[4]\t$$db{$rsname}[5]\n";}

	return($a1,$a2,$adb)
}
###################################################################################################################################################################
sub COMPARE1000G
{
   	my ($freqbool,$merge,$new,$kG,$hapmap,$dbSNPhg18_SNPrs,$freq,$path,$path_report,$path_db,$file,$Inf,$listSnpRef)=@_;
	$freqbool = shift;

	#Open all out files.
			
	my $fichier_out="$path_report/newbimfile_$Inf.out";
	my $fichier_corresp="$path_report/translation_$Inf.out";
	my $control="$path_report/control_$Inf.out";
	my $warnstrand="$path_report/warnstrand_$Inf.out";
	my $fichier_freq="$path_report/frequencies_$Inf.out";
	my $missFile="$path_report/missfile1000g_$Inf.out";
	my $outsideKG="$path_report/SNPs_outside1000g_$Inf.out";# SNPs not into 1000 genomes database
	my $exclude="$path_report/exclude_$Inf.out";

	open(NEWBIMFILE,">$fichier_out") or die "cannot open NEWBIMFILE: $fichier_out $!\n";
	open(TRANSFILE,">$fichier_corresp") or die "cannot open TRANSFILE : $fichier_corresp $!\n";
	open(WARN,">$warnstrand") or die "cannot open WARNSTRAND : $warnstrand $!\n";
	open(TROFILE,">$control") or die "cannot open TROFILE : $control $!\n";
	open(MISSFILE,">$missFile") or die "cannot open MISSFILE : $missFile $!\n";
	open(FREQFILE,">$fichier_freq") or die "cannot open FREQFILE: $fichier_freq $!\n";
	open(OUTFILE,">$outsideKG") or die "cannot open OUTFILE: $outsideKG $!\n"; #SNPs not in 1000 Genomes database
	open(EXCLUDE,">$exclude") or die "cannot open EXCLUDE: $exclude $!\n"; #SNPs to del

	#print header of out frequencies file.
	if($freqbool==1){
		print FREQFILE "chrom\trsname\tchromEnd\ta0\ta1\tfreqDB_a0\tfreqData_a0\tdelta_freq\tGenetic_Map\tblock\n"; 
	}
	else{
		print FREQFILE "chrom\trsname\tchromEnd\ta0\ta1\tfreqDB_a0\tGenetic_Map\tblock\n"; 
	}


    	open (FHI,"$file") or die "cannot open file $file : $!\n";
	
	
	my $cptSNP=0;
	my $nbNull=0;
	
	my $previousRS;
	my %previousPos=();
	my %warnMiss=();
	my $printNewbim=0;
	
  
    	while (my $line = <FHI>)
    	{
		$line =~ s/\012//g;
		$line =~ s/\015//g;
		my @elt = split(/\t/,$line);
		my $chr=$elt[0];
		my $rsname=$elt[1];
		my $posgen=$elt[2];
                my $altName=$elt[0]."-".$elt[3];
		my $pos=$elt[3];
		my $all1=$elt[4];
		my $all2=$elt[5];

		#Boolean var to allow print in Newbimfile

		my $printNewbim=0;

		#Frequency of allele observed
		my $freqDATA=$$freq{$rsname,$all1};

		my $deltafreq="";
		my $allmono;
		my $alldb;

		#Update rsname if necessary
		my $oldname=$rsname;


       		if ($$merge{$rsname}==1) 
		{
			$rsname=$$new{$rsname};
			print TRANSFILE "$oldname\tdevient\t$rsname\n";
        	}

			#Replace 0 if is a monomorphisme 
		if ($all1 eq "0")
		{
			$all1=$all2;
			$allmono=1;
		}
		
		elsif ($all2 eq "0")
		{
			$all2=$all1;
			$allmono=2;
		}
		else{}

		#Create a list of SNP to control mulplicate in 1000genomes data	
		$$listSnpRef{$rsname}=1;
		$cptSNP++;

		#Position alignement with 1000genomes datas
        		if ( $$kG{$rsname}[0]==1 )
			{	                                                                 	
				$pos=$$kG{$rsname}[3];
				$chr=$$kG{$rsname}[1];
				$posgen=$$kG{$rsname}[7] if ($$kG{$rsname}[7] ne "NA"); #By default posgen equal 0

				print TROFILE "$chr\t$rsname\t1\tNA\t$pos\n";
				
				if (($$kG{$rsname}[4] eq $all1) && ($$kG{$rsname}[5] eq $all2)) # a0db = a0data et a1db = a1data
				{
					$deltafreq=sprintf("%.4f",abs($$kG{$rsname}[6]{$all1}-$freqDATA)) if($freqbool==1);
					print FREQFILE "$chr\t$rsname\t$pos\t$all1\t$all2\t$$kG{$rsname}[6]{$all1}\t$freqDATA\t$deltafreq\t$$kG{$rsname}[7]\t$$kG{$rsname}[8]\n"; 
				}
				elsif (($$kG{$rsname}[5] eq $all1) && ($$kG{$rsname}[4] eq $all2)) # a1db = a0data et a0db = a1data
				{
					$deltafreq=sprintf("%.4f",abs($$kG{$rsname}[6]{$all2}-$freqDATA)) if($freqbool==1);
					print FREQFILE "$chr\t$rsname\t$pos\t$all1\t$all2\t$$kG{$rsname}[6]{$all2}\t$freqDATA\t$deltafreq\t$$kG{$rsname}[7]\t$$kG{$rsname}[8]\n";

				}
				elsif (($all1 eq "-") || ($all2 eq "-")) # it's a indel.
				{
					if (($all1 eq "-") && ($all2 eq "-"))
					{
						$deltafreq="NA" if($freqbool==1);
						print FREQFILE "$chr\t$rsname\t$pos\t$all1\t$all2\tNA\t$freqDATA\t$deltafreq\t$$kG{$rsname}[7]\t$$kG{$rsname}[8]\n";
					}
					elsif ($all1 eq "-")
					{		
						$deltafreq=sprintf("%.4f",abs($$kG{$rsname}[6]{$all2}-$freqDATA)) if($freqbool==1);
						print FREQFILE "$chr\t$rsname\t$pos\t$all2\t$all1\t$$kG{$rsname}[6]{$all2}\t$freqDATA\t$deltafreq\t$$kG{$rsname}[7]\t$$kG{$rsname}[8]\n";
					}
					elsif ($all2 eq "-")
					{		
						$deltafreq=sprintf("%.4f",abs($$kG{$rsname}[6]{$all1}-$freqDATA)) if($freqbool==1);
						print FREQFILE "$chr\t$rsname\t$pos\t$all1\t$all2\t$$kG{$rsname}[6]{$all1}\t$freqDATA\t$deltafreq\t$$kG{$rsname}[7]\t$$kG{$rsname}[8]\n";
					}
					

				}
				elsif (($$kG{$rsname}[4] eq STRAND($all1)) && ($$kG{$rsname}[5] eq STRAND ($all2))) # a0db = a0data(autre brin) et a1db = a1data(autre brin)
				{		
					print WARN "1\t$chr\t$rsname\t$all1\t$all2\t($$kG{$rsname}[4] $$kG{$rsname}[5]) S\n";
					my $allele_inv=STRAND($all1);
 #print "ATTENTIONS STRAND $rsname $all1 $allele_inv STRAND($all1) $$kG{$rsname}[6]{$allele_inv}\n";
					$deltafreq=sprintf("%.4f",abs($$kG{$rsname}[6]{$allele_inv}-$freqDATA)) if($freqbool==1);
					print FREQFILE "$chr\t$rsname\t$pos\t$all1\t$all2\t$$kG{$rsname}[6]{$allele_inv}\t$freqDATA\t$deltafreq\t$$kG{$rsname}[7]\t$$kG{$rsname}[8]\n";
					$all1=STRAND($all1);
					$all2=STRAND($all2);
				}
		#		elsif (($$kG{$rsname}[5] eq STRAND($all1)) && ($$kG{$rsname}[6] eq STRAND ($all2)))
				elsif (($$kG{$rsname}[5] eq STRAND($all1)) && ($$kG{$rsname}[4] eq STRAND ($all2))) # a1db = a0data(autre brin) et a0db = a1data(autre brin)
				{		
					print WARN "1\t$chr\t$rsname\t$all1\t$all2\t($$kG{$rsname}[4] $$kG{$rsname}[5]) S\n";
					my $allele_inv=STRAND($all2);
					$deltafreq=sprintf("%.4f",abs($$kG{$rsname}[6]{$allele_inv}-$freqDATA)) if($freqbool==1);
					print FREQFILE "$chr\t$rsname\t$pos\t$all1\t$all2\t$$kG{$rsname}[6]{$allele_inv}\t$freqDATA\t$deltafreq\t$$kG{$rsname}[7]\t$$kG{$rsname}[8]\n";
					$all1=STRAND($all1);
					$all2=STRAND($all2);
				}
				elsif ($all1 eq $all2)
				{	
					my($all1,$all2,$alldb)=MONOMORPH($allmono,$all1,$all2,$alldb,\%$kG,$rsname,$chr);

					$deltafreq=sprintf("%.4f",abs($$kG{$rsname}[6]{$alldb}-$freqDATA)) if($freqbool==1);
					print FREQFILE "$chr\t$rsname\t$pos\t$all1\t$all2\t$$kG{$rsname}[6]{$alldb}\t$freqDATA\t$deltafreq\t$$kG{$rsname}[7]\t$$kG{$rsname}[8]\n";
				}
				
				else
				{
					# change alternative allele to match with 1000Genomes database
					if($all1 eq $$kG{$rsname}[4]) {$alldb=$all1;$all2=$$kG{$rsname}[5];}
					elsif($all1 eq $$kG{$rsname}[5]) {$alldb=$all1;$all2=$$kG{$rsname}[4];}
					elsif($all2 eq $$kG{$rsname}[4]) {$alldb=$all2;$all1=$$kG{$rsname}[5];}
					elsif($all2 eq $$kG{$rsname}[5]) {$alldb=$all2;$all1=$$kG{$rsname}[4];}

					$deltafreq=sprintf("%.4f",abs($$kG{$rsname}[6]{$alldb}-$freqDATA)) if($freqbool==1); 
					print FREQFILE "$chr\t$rsname\t$pos\t$all1\t$all2\t$$kG{$rsname}[6]{$alldb}\t$freqDATA\t$deltafreq\t$$kG{$rsname}[7]\t$$kG{$rsname}[8]\n";
					print WARN "1\t$chr\t$rsname\t$all1\t$all2\t($$kG{$rsname}[4] $$kG{$rsname}[5]) W\n";
					
				}
				$printNewbim=1;
				
			}
			elsif ($$kG{$altName}[0]==1)
			{
			
				$pos=$$kG{$altName}[3];
				$chr=$$kG{$altName}[1];
				$posgen=$$kG{$altName}[7] if ($$kG{$altName}[7] ne "NA");
							
				print TROFILE "$chr\t$rsname\t1\t$altName\t$pos\n";
				
				if (($$kG{$altName}[4] eq $all1) && ($$kG{$altName}[5] eq $all2))
				{
					$deltafreq=sprintf("%.4f",abs($$kG{$altName}[6]{$all1}-$freqDATA)) if($freqbool==1);
					print FREQFILE "$chr\t$rsname\t$pos\t$all1\t$all2\t$$kG{$altName}[6]{$all1}\t$freqDATA\t$deltafreq\t$$kG{$altName}[7]\t$$kG{$altName}[8]\n";
				}
				elsif (($$kG{$altName}[5] eq $all1) && ($$kG{$altName}[4] eq $all2))
				{
					$deltafreq=sprintf("%.4f",abs($$kG{$altName}[6]{$all1}-$freqDATA)) if($freqbool==1);
					print FREQFILE "$chr\t$rsname\t$pos\t$all1\t$all2\t$$kG{$altName}[6]{$all1}\t$freqDATA\t$deltafreq\t$$kG{$altName}[7]\t$$kG{$altName}[8]\n";
				}
				elsif (($all1 eq "-") || ($all2 eq "-")) # The ancestral allele is unknown.
				{
					if (($all1 eq "-") && ($all2 eq "-"))
					{
						$deltafreq="NA" if($freqbool==1);
						print FREQFILE "$chr\t$rsname\t$pos\t$all1\t$all2\tNA\t$freqDATA\t$deltafreq\t$$kG{$altName}[7]\t$$kG{$altName}[8]\n";
					}
					elsif ($all1 eq "-")
					{		
						$deltafreq=sprintf("%.4f",abs($$kG{$altName}[6]{$all2}-$freqDATA)) if($freqbool==1);
						print FREQFILE "$chr\t$rsname\t$pos\t$all2\t$all1\t$$kG{$altName}[6]{$all2}\t$freqDATA\t$deltafreq\t$$kG{$altName}[7]\t$$kG{$altName}[8]\n";
					}
					elsif ($all2 eq "-")
					{	
						$deltafreq=sprintf("%.4f",abs($$kG{$altName}[6]{$all1}-$freqDATA)) if($freqbool==1);	
						if(!exists($$kG{$altName}[6]{$all1})){$$kG{$altName}[6]{$all1}="NA";$deltafreq="NA";};
						print FREQFILE "$chr\t$rsname\t$pos\t$all1\t$all2\t$$kG{$altName}[6]{$all1}\t$freqDATA\t$deltafreq\t$$kG{$altName}[7]\t$$kG{$altName}[8]\n";
					}
					
				}
				elsif (($$kG{$altName}[4] eq STRAND($all1)) && ($$kG{$altName}[5] eq STRAND ($all2)))
				{		
					print WARN "1\t$chr\t$rsname\t$all1\t$all2\t($$kG{$altName}[4] $$kG{$altName}[5]) S\n";
					my $allele_inv=STRAND($all1);
					$deltafreq=sprintf("%.4f",abs($$kG{$altName}[6]{$allele_inv}-$freqDATA)) if($freqbool==1);
					print FREQFILE "$chr\t$rsname\t$pos\t$all1\t$all2\t$$kG{$altName}[6]{$allele_inv}\t$freqDATA\t$deltafreq\t$$kG{$altName}[7]\t$$kG{$altName}[8]\n";
					$all1=STRAND($all1);
					$all2=STRAND($all2);
				}
				elsif (($$kG{$altName}[5] eq STRAND($all1)) && ($$kG{$altName}[4] eq STRAND ($all2)))
				{		
					print WARN "1\t$chr\t$rsname\t$all1\t$all2\t($$kG{$altName}[4] $$kG{$altName}[5]) S\n";
					my $allele_inv=STRAND($all2);
					$deltafreq=sprintf("%.4f",abs($$kG{$altName}[6]{$allele_inv}-$freqDATA)) if($freqbool==1);
					print FREQFILE "$chr\t$rsname\t$pos\t$all1\t$all2\t$$kG{$altName}[6]{$allele_inv}\t$freqDATA\t$deltafreq\t$$kG{$altName}[7]\t$$kG{$altName}[8]\n";
					$all1=STRAND($all1);
					$all2=STRAND($all2);
				}
				elsif ($all1 eq $all2)
				{
					my($all1,$all2,$alldb)=MONOMORPH($allmono,$all1,$all2,$alldb,\%$kG,$altName,$chr);
					$deltafreq=sprintf("%.4f",abs($$kG{$altName}[6]{$all1}-$freqDATA)) if($freqbool==1);
					print FREQFILE "$chr\t$rsname\t$pos\t$all1\t$all2\t$$kG{$altName}[6]{$all1}\t$freqDATA\t$deltafreq\t$$kG{$altName}[7]\t$$kG{$altName}[8]\n";
				}
				else
				{	
					my $alldb;

					# change alternative allele to match with 1000Genomes database
					if($all1 eq $$kG{$altName}[4]) {$alldb=$all1;$all2=$$kG{$altName}[5];}
					elsif($all1 eq $$kG{$altName}[5]) {$alldb=$all1;$all2=$$kG{$altName}[4];}
					elsif($all2 eq $$kG{$altName}[4]) {$alldb=$all2;$all1=$$kG{$altName}[5];}
					elsif($all2 eq $$kG{$altName}[5]) {$alldb=$all2;$all1=$$kG{$altName}[4];}

					$deltafreq=sprintf("%.4f",abs($$kG{$altName}[6]{$alldb}-$freqDATA)) if($freqbool==1);
					print FREQFILE "$chr\t$rsname\t$pos\t$all1\t$all2\t$$kG{$altName}[6]{$alldb}\t$freqDATA\t$deltafreq\t$$kG{$altName}[7]\t$$kG{$altName}[8]\n";
					print WARN "1\t$chr\t$rsname\t$all1\t$all2\t($$kG{$altName}[4] $$kG{$altName}[5]) W\n";
				}
				
				$printNewbim=1;
			
			}
			elsif ( $$kG{$oldname}[0]==1 )
			{	                                          
    	                                         	
				$pos=$$kG{$oldname}[3];
				$chr=$$kG{$oldname}[1];
				$posgen=$$kG{$oldname}[7] if ($$kG{$oldname}[7] ne "NA");

				print TROFILE "$chr\t$oldname\t1\tNA\t$pos\n";

				if (($$kG{$oldname}[4] eq $all1) && ($$kG{$oldname}[5] eq $all2))
				{
					$deltafreq=sprintf("%.4f",abs($$kG{$oldname}[6]{$all1}-$freqDATA)) if($freqbool==1);
					print FREQFILE "$chr\t$oldname\t$pos\t$all1\t$all2\t$$kG{$oldname}[6]{$all1}\t$freqDATA\t$deltafreq\t$$kG{$oldname}[7]\t$$kG{$oldname}[8]\n";
				}
				elsif (($$kG{$oldname}[5] eq $all1) && ($$kG{$oldname}[4] eq $all2))
				{
					$deltafreq=sprintf("%.4f",abs($$kG{$oldname}[6]{$all1}-$freqDATA)) if($freqbool==1);
					print FREQFILE "$chr\t$oldname\t$pos\t$all1\t$all2\t$$kG{$oldname}[6]{$all1}\t$freqDATA\t$deltafreq\t$$kG{$oldname}[7]\t$$kG{$oldname}[8]\n";
				}
				elsif (($all1 eq "-") || ($all2 eq "-")) # it's a indel.
				{
					if (($all1 eq "-") && ($all2 eq "-"))
					{
						$deltafreq="NA" if($freqbool==1);
						print FREQFILE "$chr\t$oldname\t$pos\t$all1\t$all2\tNA\t$freqDATA\t$deltafreq\t$$kG{$oldname}[7]\t$$kG{$oldname}[8]\n";
					}					
					elsif ($all1 eq "-")
					{		
						$deltafreq=sprintf("%.4f",abs($$kG{$oldname}[6]{$all2}-$freqDATA)) if($freqbool==1);
						print FREQFILE "$chr\t$oldname\t$pos\t$all2\t$all1\t$$kG{$oldname}[6]{$all2}\t$freqDATA\t$deltafreq\t$$kG{$oldname}[7]\t$$kG{$oldname}[8]\n";
					}
					elsif ($all2 eq "-")
					{	
						$deltafreq=sprintf("%.4f",abs($$kG{$oldname}[6]{$all1}-$freqDATA)) if($freqbool==1);
						print FREQFILE "$chr\t$oldname\t$pos\t$all1\t$all2\t$$kG{$oldname}[6]{$all1}\t$freqDATA\t$deltafreq\t$$kG{$oldname}[7]\t$$kG{$oldname}[8]\n";
					}
					
				}
				elsif (($$kG{$oldname}[4] eq STRAND($all1)) && ($$kG{$oldname}[5] eq STRAND ($all2)))
				{		
					print WARN "1\t$chr\t$oldname\t$all1\t$all2\t($$kG{$oldname}[4] $$kG{$oldname}[5]) S\n";
					my $allele_inv=STRAND($all1);
					$deltafreq=sprintf("%.4f",abs($$kG{$oldname}[6]{$allele_inv}-$freqDATA)) if($freqbool==1);
					print FREQFILE "$chr\t$oldname\t$pos\t$all1\t$all2\t$$kG{$oldname}[6]{$allele_inv}\t$freqDATA\t$deltafreq\t$$kG{$oldname}[7]\t$$kG{$oldname}[8]\n";
					$all1=STRAND($all1);
					$all2=STRAND($all2);
				}
				elsif (($$kG{$oldname}[5] eq STRAND($all1)) && ($$kG{$oldname}[4] eq STRAND ($all2)))
				{		
					print WARN "1\t$chr\t$oldname\t$all1\t$all2\t($$kG{$oldname}[4] $$kG{$oldname}[5]) S\n";
					my $allele_inv=STRAND($all2);
					$deltafreq=sprintf("%.4f",abs($$kG{$oldname}[6]{$allele_inv}-$freqDATA)) if($freqbool==1);
					print FREQFILE "$chr\t$oldname\t$pos\t$all1\t$all2\t$$kG{$oldname}[6]{$allele_inv}\t$freqDATA\t$deltafreq\t$$kG{$oldname}[7]\t$$kG{$oldname}[8]\n";
					$all1=STRAND($all1);
					$all2=STRAND($all2);
				}
				elsif ($all1 eq $all2)
				{
					my($all1,$all2,$alldb)=MONOMORPH($allmono,$all1,$all2,$alldb,\%$kG,$oldname,$chr);
					$deltafreq=sprintf("%.4f",abs($$kG{$oldname}[6]{$all1}-$freqDATA)) if($freqbool==1);
					print FREQFILE "$chr\t$oldname\t$pos\t$all1\t$all2\t$$kG{$oldname}[6]{$all1}\t$freqDATA\t$deltafreq\t$$kG{$oldname}[7]\t$$kG{$oldname}[8]\n";
				}
				else
				{
					my $alldb;
		
					# change allele to match with 1000Genomes database
					if($all1 eq $$kG{$oldname}[4]) {$alldb=$all1;$all2=$$kG{$oldname}[5];}
					elsif($all1 eq $$kG{$oldname}[5]) {$alldb=$all1;$all2=$$kG{$oldname}[4];}
					elsif($all2 eq $$kG{$oldname}[4]) {$alldb=$all2;$all1=$$kG{$oldname}[5];}
					elsif($all2 eq $$kG{$oldname}[5]) {$alldb=$all2;$all1=$$kG{$oldname}[4];}

					$deltafreq=sprintf("%.4f",abs($$kG{$oldname}[6]{$alldb}-$freqDATA)) if($freqbool==1);
					print FREQFILE "$chr\t$oldname\t$pos\t$all1\t$all2\t$$kG{$oldname}[6]{$alldb}\t$freqDATA\t$deltafreq\t$$kG{$oldname}[7]\t$$kG{$oldname}[8]\n";
					print WARN "1\t$chr\t$oldname\t$all1\t$all2\t($$kG{$oldname}[4] $$kG{$oldname}[5]) W\n";
					
				}
				$printNewbim=1;	
				
			}
			elsif ($$hapmap{$rsname}[0]==1)
			{
			 
				$pos=$$hapmap{$rsname}[3];
				$chr=$$hapmap{$rsname}[1];
				$posgen=$$hapmap{$rsname}[7] if ($$hapmap{$rsname}[7] ne "NA");

				print TROFILE "$chr\t$rsname\t2\tNA\t$pos\n";
				print OUTFILE "chr$chr\t$pos\t".($pos+1)."\t$rsname\n";
				if (($$hapmap{$rsname}[4] eq $all1) && ($$hapmap{$rsname}[5] eq $all2))
				{
					$deltafreq=sprintf("%.4f",abs($$hapmap{$rsname}[6]{$all1}-$freqDATA)) if($freqbool==1);
					print FREQFILE "$chr\t$rsname\t$pos\t$all1\t$all2\t$$hapmap{$rsname}[6]{$all1}\t$freqDATA\t$deltafreq\t$$hapmap{$rsname}[7]\t$$hapmap{$rsname}[8]\n";
				}
				elsif (($$hapmap{$rsname}[5] eq $all1) && ($$hapmap{$rsname}[4] eq $all2))
				{
					$deltafreq=sprintf("%.4f",abs($$hapmap{$rsname}[6]{$all1}-$freqDATA)) if($freqbool==1);
					print FREQFILE "$chr\t$rsname\t$pos\t$all1\t$all2\t$$hapmap{$rsname}[6]{$all1}\t$freqDATA\t$deltafreq\t$$hapmap{$rsname}[7]\t$$hapmap{$rsname}[8]\n";
				}
				elsif (($all1 eq "-") || ($all2 eq "-")) # it's a indel.
				{
					if ( ($all1 eq "-") && ($all2 eq "-"))
					{
						$deltafreq="NA" if($freqbool==1);
						print FREQFILE "$chr\t$rsname\t$pos\t$all1\t$all2\tNA\t$freqDATA\t$deltafreq\t$$hapmap{$rsname}[7]\t$$hapmap{$rsname}[8]\n";
					}
					
					elsif ($all1 eq "-")
					{		
						$deltafreq=sprintf("%.4f",abs($$hapmap{$rsname}[6]{$all2}-$freqDATA)) if($freqbool==1);
						print FREQFILE "$chr\t$rsname\t$pos\t$all1\t$all2\t$$hapmap{$rsname}[6]{$all2}\t$freqDATA\t$deltafreq\t$$hapmap{$rsname}[7]\t$$hapmap{$rsname}[8]\n";
					}
					elsif ($all2 eq "-")
					{		
						$deltafreq=abs($$hapmap{$rsname}[6]{$all1}-$freqDATA) if($freqbool==1);
						print FREQFILE "$chr\t$rsname\t$pos\t$all1\t$all2\t$$hapmap{$rsname}[6]{$all1}\t$freqDATA\t$deltafreq\t$$hapmap{$rsname}[7]\t$$hapmap{$rsname}[8]\n";
					}
					
				}
				elsif (($$hapmap{$rsname}[4] eq STRAND($all1)) && ($$hapmap{$rsname}[5] eq STRAND ($all2)))
				{		
					print WARN "1\t$chr\t$rsname\t$all1\t$all2\t($$hapmap{$rsname}[4] $$hapmap{$rsname}[5]) S\n";
					my $allele_inv=STRAND($all1);
					$deltafreq=sprintf("%.4f",abs($$hapmap{$rsname}[6]{$allele_inv}-$freqDATA)) if($freqbool==1);
					print FREQFILE "$chr\t$rsname\t$pos\t$all1\t$all2\t$$hapmap{$rsname}[6]{$allele_inv}\t$freqDATA\t$deltafreq\t$$hapmap{$rsname}[7]\t$$hapmap{$rsname}[8]\n";
					$all1=STRAND($all1);
					$all2=STRAND($all2);
				}
				elsif (($$hapmap{$rsname}[5] eq STRAND($all1)) && ($$hapmap{$rsname}[4] eq STRAND ($all2)))
				{		
					print WARN "1\t$chr\t$rsname\t$all1\t$all2\t($$hapmap{$rsname}[4] $$hapmap{$rsname}[5]) S\n";
					my $allele_inv=STRAND($all2);
					$deltafreq=sprintf("%.4f",abs($$hapmap{$rsname}[6]{$allele_inv}-$freqDATA)) if($freqbool==1);
					print FREQFILE "$chr\t$rsname\t$pos\t$all1\t$all2\t$$hapmap{$rsname}[6]{$allele_inv}\t$freqDATA\t$deltafreq\t$$hapmap{$rsname}[7]\t$$hapmap{$rsname}[8]\n";
					$all1=STRAND($all1);
					$all2=STRAND($all2);
				}
				elsif ($all1 eq $all2)
				{
					my($all1,$all2,$alldb)=MONOMORPH($allmono,$all1,$all2,$alldb,\%$hapmap,$rsname,$chr);
					$deltafreq=sprintf("%.4f",abs($$hapmap{$rsname}[6]{$all1}-$freqDATA)) if($freqbool==1);
					print FREQFILE "$chr\t$rsname\t$pos\t$all1\t$all2\t$$hapmap{$rsname}[6]{$all1}\t$freqDATA\t$deltafreq\t$$hapmap{$rsname}[7]\t$$hapmap{$rsname}[8]\n";
				}
				else
				{
					my $alldb;
					# change alternative allele to match with hapmap database
					if($all1 eq $$hapmap{$rsname}[4]) {$alldb=$all1;$all2=$$hapmap{$rsname}[5];}
					elsif($all1 eq $$hapmap{$rsname}[5]) {$alldb=$all1;$all2=$$hapmap{$rsname}[4];}
					elsif($all2 eq $$hapmap{$rsname}[4]) {$alldb=$all2;$all1=$$hapmap{$rsname}[5];}
					elsif($all2 eq $$hapmap{$rsname}[5]) {$alldb=$all2;$all1=$$hapmap{$rsname}[4];}

					$deltafreq=sprintf("%.4f",abs($$hapmap{$rsname}[6]{$alldb}-$freqDATA)) if($freqbool==1);
					print FREQFILE "$chr\t$rsname\t$pos\t$all1\t$all2\t$$hapmap{$rsname}[6]{$alldb}\t$freqDATA\t$deltafreq\t$$hapmap{$rsname}[7]\t$$hapmap{$rsname}[8]\n";
					print WARN "1\t$chr\t$rsname\t$all1\t$all2\t($$hapmap{$rsname}[4] $$hapmap{$rsname}[5]) W\n";
					
				}
				
				$printNewbim=1;	 
					
			}		
			else    # then we see on dbSNP 
			{
				
				$dbSNPhg18_SNPrs->execute($rsname);
				if ($dbSNPhg18_SNPrs->rows!=0)
				{
					my @result_query=$dbSNPhg18_SNPrs->fetchrow_array;
					$chr=$result_query[0];
					$chr=~s/chrX/23/g;
					$chr=~s/chrY/24/g;
					$chr=~s/chrM/26/g;
					$chr=~s/chr(\d+).*/\1/g;
					$pos=$result_query[1];
					print TROFILE "$chr\t$rsname\t3\tNA\t$pos\n";
					print OUTFILE "chr$chr\t$pos\t".($pos+1)."\t$rsname\n";
					$deltafreq="NA" if($freqbool==1);
					print FREQFILE "$chr\t$rsname\t$pos\t$all1\t$all2\tNA\t$freqDATA\t$deltafreq\tNA\tNA\n";
					$printNewbim=1;	
				
				}		
				else
				{
					  $warnMiss{$rsname}[0]=1; #we flag the SNP ligne isn't match in database
					  $warnMiss{$rsname}[1]=$chr;
					  $warnMiss{$rsname}[2]=$rsname;
					  $warnMiss{$rsname}[3]=$previousPos{$previousRS}; #previous position
					  $warnMiss{$rsname}[4]=$all1;
					  $warnMiss{$rsname}[5]=$all2;
					  $warnMiss{$rsname}[6]=$freqDATA; #frequency of allele A1
						
					  $nbNull++;
					  print MISSFILE "$chr\t$rsname\n";
						
				}
				
			}
				

			if ( $warnMiss{$previousRS}[0]==1 )
			{
				my $newpos=int(($pos+$warnMiss{$previousRS}[3])/2);
				print TROFILE "$warnMiss{$previousRS}[1]\t$warnMiss{$previousRS}[2]\t4\tNA\t$newpos\n";
				print OUTFILE "chr$warnMiss{$previousRS}[1]\t$newpos\t".($newpos+1)."\t$warnMiss{$previousRS}[2]\n";
				print FREQFILE "$warnMiss{$previousRS}[1]\t$warnMiss{$previousRS}[2]\t$newpos\t$warnMiss{$previousRS}[4]\t$warnMiss{$previousRS}[5]\tNA\t$warnMiss{$previousRS}[6]\tNA\tNA\tNA\n";
				print NEWBIMFILE "$warnMiss{$previousRS}[1]\t$warnMiss{$previousRS}[2]\t0\t$newpos\t$warnMiss{$previousRS}[4]\t$warnMiss{$previousRS}[5]\n";

			}   	                                

			print NEWBIMFILE "$chr\t$rsname\t$posgen\t$pos\t$all1\t$all2\n" if( $printNewbim==1 );    	        	
			$previousPos{$rsname}=$pos;
			$previousRS=$rsname;
			
	}

		close(NEWBIMFILE);
		close(TRANSFILE);
		close(TROFILE);
		close(WARNSTRAND);
		close(FREQFILE);
		close(MISSFILE);
		close(OUTFILE);
		close(EXCLUDE);

		close (FHI);
		return ($cptSNP,$nbNull);
		
 
}
###################################################################################################################################################################

sub CONTROL_BUILD
{
	my ($datafile,$dbSNPhg17_SNPpos,$dbSNPhg18_SNPpos,$dbSNPhg19_SNPpos)=@_;
	my %file=();
	my %nbdiff=();
        my $i=0;
	my $index;
	
	my $error17=0;
	my $error18=0;
	my $error19=0;
	
	
	open (FILE,"$datafile") or die " cannot open Data file (CONTROLE BUILD) :$!\n";
	my @lines=<FILE>;
	srand;

	my $index;
	for(my $i=0;$i<500;$i++)
	{
		$index=rand($#lines);
		chomp ($lines[$index]);
		my($chrom,$rs,$pos,$chromEnd,$a,$b)=split(/\t/,$lines[$index]);
		$file{$rs}=$chromEnd;
		splice (@lines, $index, 1);
		
	}
	@lines=();
	close(FILE);
	
		

		foreach my $rs (keys %file)
		{ 
  			$dbSNPhg17_SNPpos->execute($rs);
			my $pos17 = $dbSNPhg17_SNPpos->fetchrow_array;
  			if ($pos17!=$file{$rs}){$error17++}; 
			$dbSNPhg18_SNPpos ->execute($rs);
  			my $pos18 = $dbSNPhg18_SNPpos->fetchrow_array;
  			if ($pos18!=$file{$rs}){$error18++}; 
			$dbSNPhg19_SNPpos->execute($rs);
  			my $pos19 = $dbSNPhg19_SNPpos->fetchrow_array;
  			if ($pos19!=$file{$rs}){$error19++}; 
  
		}
  	
	if ($error19<50){ return 19;}
	elsif  ($error18<50){return 18;}
	elsif  ($error17<50){ return 17;}
	else{ print "\n*******Build unknown*******\n"; return 0}
		
}

###################################################################################################################################################################
sub STRAND
{ 
	   my ($allele)=@_;
	   if ($allele eq "A") { return "T"; }
	elsif ($allele eq "T") { return "A"; }
	elsif ($allele eq "C") { return "G"; }
	elsif ($allele eq "G") { return "C"; }

}
###################################################################################################################################################################
sub UPDATERSNAME
{
	my ($PATH,$dbRef,$buildref,$pathdata,$pathfreq,$filename,$out,$date)= @_;
   
	print "\n*******PROCESSING UPDATE RSNAME FOR FILE: $filename *******\n";
   	
	my $dbDIR=$$dbRef{"path"};

	my $MergeFile = "$dbDIR/RsMergeArch.bcp";
	my %existMerge=();
	my %newName=();

	# the file of frequencies has a name like bim file but with suffix .frq ( we match the good file of frequencies )  
	my $freqfile=$filename;
	$freqfile=~ s/(.*)\.bim/\1\.frq/g;


	#my $tmpRep="$PATH/tmp"; # Make tempory programme's workspace. It's emptied when the programme was terminated
  	#MAKEDIRE($tmpRep);
   
#                           /////////////////////////////////////////////////////////////////////////////////////////
#                                                  LOAD ON MEMORY RsMERGEARCH FILE 
  	print "1-LOAD RSMERGE DATAFILE\n";
	LOAD_RsMERGE($MergeFile,\%newName,\%existMerge); 
	my $freqb=0;#boolean 

#				/////////////////////////////////////////////////////////////////////////////////////////
#                                                 LOAD FREQUENCIES DATA  
	my %freqDATA=();
	my $data="$pathfreq/$freqfile";
	if ( -e "$pathfreq/$freqfile" ) #it's optional 
	{
		$freqb=1;
		print "2-LOAD FREQUENCY DATAFILE : $data\n";
		FREQ_DATA($data,\%freqDATA);
	}
#                           /////////////////////////////////////////////////////////////////////////////////////////
#                                                  UPDATE BLOCK DEFINITION WITH RECOMBINATION RATE THRESHOLD DEFINED BY USER  

	
	
	my $path1kG=$$dbRef{"Path1kg"};
	my $pathHapmap=$$dbRef{"PathHapmap"};
	my $pre=$$dbRef{"RecMap_PRE"};
	my $post1=$$dbRef{"RecMap_POST1"};
	my $post2=$$dbRef{"RecMap_POST2"};
	my $RecRate=$$dbRef{"Recombination_Rate"};
	
	
	my $blockDataOut1kg=$path1kG."/".$pre."1000G".$post1.".annot";
	print "3-UPDATE BLOCK DEFINITION FILES ON 1000G WITH RECOMBINATION RATE THRESHOLD DEFINED BY USER : $RecRate\n";
	UPDATE_RM($path1kG,$pre,$post1,$RecRate,$blockDataOut1kg);

	my $blockDataOutHapmap=$pathHapmap."/".$pre."Hapmap".$post2.".annot";
	print "4-UPDATE BLOCK DEFINITION FILES ON HAPMAP WITH RECOMBINATION RATE THRESHOLD DEFINED BY USER : $RecRate\n";
	UPDATE_RM($pathHapmap,$pre,$post2,$RecRate,$blockDataOutHapmap);




#                           /////////////////////////////////////////////////////////////////////////////////////////
#                                                  LOAD ON MEMORY 1000GENOMES AND HAPMAP PSEUDO-DATABASE  
	my %kG=();
	my %hapmap=();
	
	my $KG=$$dbRef{"_1kg"};
	my $Hapmap=$$dbRef{"Hapmap"};

	print "5-LOAD 1000G DATAFILE\n";
	LOAD_DB(\%kG,$dbDIR ,$blockDataOut1kg,$KG);
	print "6-LOAD HAPMAP DATAFILE\n";
	LOAD_DB(\%hapmap, $dbDIR,$blockDataOutHapmap,$Hapmap);
	

#                           /////////////////////////////////////////////////////////////////////////////////////////
#					CONNECT AND PREPARE STATEMENTS ON LOCAL DBSNP DATABASE 

	my $dsn="hg19";
	my $host="172.18.254.161";
	my $user="anonymous"; #"aduval";
        my $password=""; ###"aduval1984";
	# Connect to the database.
        my $dbh= DBI->connect("DBI:mysql:host=$host;database=$dsn",$user, $password,{'RaiseError' => 1});
	

	my $dbSNPhg17_SNPpos = $dbh->prepare("select chromEnd from hg17.snp125 WHERE name=?");
	my $dbSNPhg18_SNPpos = $dbh->prepare("select chromEnd from hg18.snp130 WHERE name=?");
	my $dbSNPhg19_SNPpos = $dbh->prepare("select chromEnd from hg19.snp131 WHERE name=?");
	my $dbSNPhg18_SNPrs = $dbh->prepare("select chrom,chromEnd from hg18.snp130 WHERE name=?");
		
	
#                           /////////////////////////////////////////////////////////////////////////////////////////
	
	
	my $file ="$pathdata/$filename";
	my $information =$filename."_".$date;

		print "7-CONTROLE BUILD FILE : $file\n";
		my $build=CONTROL_BUILD($file,$dbSNPhg17_SNPpos,$dbSNPhg18_SNPpos,$dbSNPhg19_SNPpos);
		print ">BUILD : hg$build\n";
	
		my %listSNP;
		my $cptSNP;

		if ($build>=17)
		{
			print "8-START COMPARE 1000G\n";
			my ($cptSNP,$nbNull)=COMPARE1000G($freqb,\%existMerge,\%newName,\%kG,\%hapmap,$dbSNPhg18_SNPrs,\%freqDATA,$PATH,$out,$dbDIR,$file,$information,\%listSNP,\$cptSNP);
			print "9-CONTROL DUPLICATE 1000G\n";
			my $nbMult=CONTROL_DUPLICATE_1000G(\%listSNP,$dbDIR,$KG,$out,$information);
			print "\n*******************************INFORMATIONS******************************\n";
			print ">NBRE SNPS $cptSNP\n";
			print ">$nbNull SNPs are not in database (build hg18) et $nbMult are more 1 copy\n";
			
		}
		else
		{
			print "\n*******PROCESS FAILED : Build hg$build not compatible *******\n";
		}

#Close all statements

		$dbSNPhg17_SNPpos->finish();
		$dbSNPhg18_SNPpos->finish();
		$dbSNPhg18_SNPrs->finish();
		$dbSNPhg19_SNPpos->finish();

# Disconnect from the local dbSNP database.
		

  		$dbh->disconnect();
return;

} #end sub
################################################################################################################################################################### 

sub CONTROL_DUPLICATE_1000G
{
	#This subroutine won't be necessary because in 1000genome release datas of june2010 the duplicate and multiplicate are deleted. Just for Control !!! 
	
	my ($listMarker,$dbDIR,$dbfile,$outdir,$inf)=@_; #$listMarker is newbimfile.out

	my $KG="$dbDIR/$dbfile";
	my%cpt=();
	my $cptSNP=0;
	my  $nbMult=0;
	my $nbNull=0;

	my $doubleFile="$outdir/double1000g_$inf.out";
	
	open (FHI, "$KG") or die "cannot open file $KG";
	open(DOUFILE,">$doubleFile") or die "Cannot open DOUBLE 1000G $doubleFile : $!\n" ;
	
	while (my $line = <FHI>) 
	{
		$line=~ s/\012//g;
		$line=~ s/\015//g;
		$line=~ s/\s+/\t/g;

		my ($chrom,$rs,$pos,$a,$b,$freq) = split(/\t/,$line); 

			if ($$listMarker{$rs}==1) 
			{
			
				if (!exists($cpt{$rs}))  #it's first time we've seen this one
				{
					$cpt{$rs} = 1 ;
    				}
				elsif ($cpt{$rs}==1) 
				{
        				$cpt{$rs}++;
    				} 
				elsif($cpt{$rs}>1)
				{      
       					$cpt{$rs}++;
        				print DOUFILE "$chrom $rs $cpt{$rs}\n";
					$nbMult++;
        			}
				else{}
			}
			
			
	}
	close(FHI);
	close(DOUFILE);
 
	return $nbMult;
	
}
###################################################################################################################################################################
1;
